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Abstract 

This is the first in a series of papers which ultimately aims on improving on 
the present estimates on the axion mass by modelling the topological non- 
perturbative QCD dynamics. Axions couple to instantons and their mass is 
set by the topological susceptibility whose temperature dependence we estimate 
with the interacting instanton liquid model (IILM). Since accurate finite tem- 
perature instanton calculations have problems and do not consider fluctuations 
in the topological charge, we develop an improved grand canonical version of 
the IILM to study topological fluctuations in the quark gluon plasma. In this 
first paper we will calibrate the model against the topological susceptibility at 
zero temperature, in the chiral regime of physical quark masses. 



1. Introduction 

The strong interactions at finite temperature are believed to display a num- 
ber of interesting, non-perturbative phenomena, among which are the confine- 
ment /deconfinement transition, spontaneous P and CP violation and chiral 
symmetry restoration. The latter is believed to have its origin in topolog- 
ical fluctuations. Lattice simulations, e.g. [H, IHij ]. and phenomenology, e.g. 
[HL EH SE Ell , have shown that the chiral dynamics of QCD is well described 
by instanton models. 

Another interesting question related to topological fluctuations is the de- 
termination of the axion mass. Axions couple to instantons and their mass is 
directly proportional to the topological susceptibility. The latter turns out to 
be a chiral property of QCD and can thus also be expected to be well described 
by the IILM. The main physical question that we want to address is the compu- 
tation of the topological susceptibility and the systematic effects that pertain to 
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its determination with instanton methods in the chiral regime of light, physical 
quark masses. 

Instanton models are based on a combination of semi-classical expansion 



47j and variational approach [13j, |14[ . Taking this variational path integral as 
a starting point, Shuryak investigated what has become known as the interact- 
ing instanton liquid model (IILM) [43, 44 1. In [3j| many bulk properties were 



computed and seen to be consistent with the available lattice data and phe- 
nomenology. Some recent studies [2(| 0] corroborate the earlier results that the 
IILM rather accurately describes the chiral properties of QCD, i.e. that instan- 
tons are the dominant degrees of freedom as far as the chiral regime of QCD is 
concerned. However, the IILM fails to reproduce confinement. 

The topological susceptibility is a key parameter of QCD and has been in- 
vestigated in many lattice studies. Comparatively few studies have addressed 
this quantity within the IILM [46j]. One reason is that, so far, the IILM is based 
on a canonical ensemble. Although one can extract the topological susceptibil- 
ity from the canonical ensemble through the decay of correlators, e.g. see [46| 
for IILM and [JJ @] for lattice simulations, it is most natural to use the grand 
canonical ensemble to study the topological susceptibility. Recent investigations 
of the IILM in the grand canonical ensemble O 53] were based on canonical 
simulations and a fugacity expansion, while we will set up a grand canonical 
IILM that uses grand canonical Monte Carlo simulatio ns: see also [lO] for a 



'mean-field' study of the grand-canonical ensemble and |35j for an exploratory 
investigation of grand canonical simulations in a simpler framework. 

While developing the gran d canonical IILM, we have found that the existing 
finite temperature IILM, [451 ] . displays an unphysical behaviour in that it does 
not allow for a thermodynamic limit. Specifically the instanton-instanton inter- 



action, Eq. (3.11) in [45[, contains a term that decays very slowly with instanton 
separation R, 

and is not integrable. In their original paper the authors do discuss this long- 
range interaction and point out that they found the 0(1/ R) dyon-dyon be- 
haviour for a wide range of intermediate separations. It might well be that 
the interactions are still well described by this ansatz for the simulation boxes 
used in subsequent numerical investigations, e.g. [39j, but for studying the large 
volume behaviour it is not appropriate. 

We remedy this deficiency by re-deriving the interactions and setting up 
a numerical framework that avoids using parametric fitting formulas, such as 
|T]); instead we will integrate the Lagrangian density for a pair of instantons 
exactly, i.e. numerically. To the best of our knowledge, the exact action density 
for a pair has not been published in the literature before; we will provide it 
for Harrington-Shepard calorons 25]. The explicit form allows us to perform 
the numerical integration in an efficient way by exploiting the symmetries of 
the integrand, and an exact analytic computation for widely separated pairs 
is possible because of the localised nature of the integrands. Hence, the large 



2 



separation interactions are under very good control. In particular, the large 
separation instanton-instanton interaction at finite temperature is not given 



by {TJ), as we will see more explicitly in the final paper of this series 54j ; 
this paper, however, we will restrict ourselves to zero temperature. We hope 
that this framework can be build upon to include the non-trivial holonomy 
calorons [2^, [3(1 (HI, which may play part in the confinement /deconfinement 
phase transition, and for which good fitting formulas will be even harder to 
come by because of their more complicated structure. As mentioned before, in 
this series we will restrict ourselves to the Harrington-Shepard calorons. 

In section [5] we will review the standard strategy used to derive the parti- 
tion function for an ensemble of background gauge fields, i.e. the semi-classical 
approximation. We will then re-derive the interactions for the so-called ratio 
ansatz, used to construct multi-instanton backgrounds from individual instan- 
tons, in section [3] and compare it with other available ansatze. In section |4] we 
present the numerical framework we have set up to deal with the simulations. 
Given that different ansatze are available, we will study their effect on some 
bulk properties in section [S] and we endeavour to get a handle on systematic un- 
certainties inherent in this approach. Finally we fix the free parameters of the 
model and summarise our results in section [SJ Finite temperature simulations 



will be dealt with in [53j and [54 1 



2. Saturating the path integral 

The IILM path integral is an approximation to the fundamental path integral 
by saturating the latter with a given ansatz for the multi-instanton background. 
The functional measure consists of small fluctuations around that classical con- 
figuration. To make analytical progress, the action is expanded to quadratic 
order to define the 'free' part that is used in perturbation theory. In general 
the background induces non-Gaussian fluctuations that need to be treated ex- 
actly. The directions of these zero modes can (sometimes) be integrated up, 
and correspond to the tangent space of a generically non-trivial manifold. The 
coordinates on this so-called moduli-space can be interpreted as those degrees 
of freedom whose quantum mechanics approximates the low-energy dynamics 
of the fundamental theory. 

In order to discuss the approximations that are eventually used, we will 
now sketch the construction of the variational path integral, paying particular 
attention to the low lying modes. Details pertaining to the variational approach, 
gauge fixing and renormalisation can be found in the original papers HI, 14 1. 



We denote by <j> the collection of bosonic fields. The classical action we write as 
S c and the classical interaction is defined as 

Sint =S C - NS , (2) 

where N is the number of instanton constituents of the background and So the 
action of an individual instanton. Assume that the background has N y (quasi) 
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zero modes, which we denote collectively by 7. We can then writ^E using the 
eigenfunctions of the free part of the action S 2 S/5<j) 2 \ , _ , at zero and finite 7, 



00 

0c(x,O) + ^2( n f] n (x,0) , 
n=l 

OO 

c (x, 7 ) + (nVn(x,l) + 0(7 2 ) . 

n=AT 7 + l 



(3) 



This can be rearranged to (omitting the x-dependence for notational clarity) 



--N-, 



OO 

E 



(TX Ul) ~ 0c(O) + J2 CnVnil) rj m (0) 



(4) 



where we used the fact that r](x, 0) forms a complete basis. Clearly ^(7 = 0, J_ 
?ji with i = 1,... ,Ny. The Jacobian for the variable change {£„} ~> {7)Cm}> 



follows from the following partial derivatives 

3mn • 



d-fi 

dCn 



7=0 



(5) 
(6) 



7=0 



Note the occurrence of the <fi part. This will lead to new interactions which have 
no classical counterpart but are purely quantum mechanical; to 1-loop order, 
we are allowed to discard them. From this matrix structure it follows that the 
Jacobian is given by det(J d ln (pr] m ). 

Now, we do not know the set {rj} of exact low lying eigenfunctions. However, 
we can approximate it by constructing an orthonormal set of the known single 
particle zero modes that descend from the exact solutions used to build up the 
background field. With a slight abuse of notation, we substitute 77 — > fj = ObV> 
Ob is the matrix that generates an orthonormal basis from the original set {rf\ 
of single particle zero modes. The Jacobian is then given by 



det 



d~. 



> Vn 



= det 



d Jn <$>f]m det B 



(7) 



It corresponds to the quantum mechanical gluonic interactions. The high- 
frequency eigenvalues, encoded in the determinant of the fluctuation operator 
S 2 S/S(j) 2 are assumed to be AT-fold degenerate, and so the fluctuation determi- 
nant factorises. 



lr This follows [13. 
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In QCD we also need to introduce quarks and treat their interactions with 
the background field. In the case where the Dirac operator admits quasi-zero 
modes we can approximate the low frequency part in the same way as for the 
gluonic case; the high-frequency fluctuations will again be assumed to factorise. 
As for the Jacobian, we do not know the exact set of low-lying eigenfunctions 
for the superposition, but we approximate it by constructing an orthonormal 
set of the exact single particle zero modes £ n , i.e. £ = Of£- The Dirac operator, 
truncated to that subspace, is then given by 

(p + m) low = W+™jij 16) ® (Ci I - (Pij + mSij) |6) ® (£■ | , (8) 

with p t j — (£i\p\£j)- The matrix of overlaps is related to the single particle 
zero mode overlaps by 

pi ow = P = OlpO F , (9) 

with pij = (£i\P\£j)- Note that this is not a similarity transformation because 
Of is not unitary. 



3. Interactions in the IILM 



We will now turn to instantons in QCD. In this paper, we will only discuss 
BPST instantons 2] . In terms of the 't Hooft potential 

n(x,{y,p}) = ^, (10) 
with r 2 — (x — y) 2 , the BPST instanton in singular gauge is given by 

\a s~\ab /-b 

d„U(x,{y,p}) 
1 + n{x,{y lP }) 

with C^ v = ffi for instantons, C„ = rf^ for anti-instantons and rj the 't Hooft 
symbols. The collective coordinates are: y the centre, p the size and O the 
colour orientation in the adjoint representation. 

The simplest background configuration is the sum ansatz, as used for in- 



stance in [13| . It was shown in 41 1, that the sum ansatz produces an unphysical 
amount of repulsion; this is due to the fact that the field strength actually 
diverges at the individual centres, and is in sharp contrast to the individual sin- 
gular gauge instanton whose field strength is finite at the central- In this case, 



2 Actually 1 + il is the 't Hooft potential. 

3 Note, however, that the field strength of the individual singular gauge instanton is not 
continuous at the centre and is only defined on the punctured Euclidean space. Incidentally, 
the winding about this singular point corresponds to the winding at infinity of the regular 
instanton. 
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Ansatz 



Description 



Re 



Interactions for ratio ansatz as derived in this paper. 



Rh 



Gluonic interactions are derived from the ratio ansatz whereas 
the quark overlaps use a sum ansatz, [39j |. 



S 



Interactions have been derived from the so-called streamline 
ansatz. These are only available at zero temperature, [39( |48| . 



Table 1: Several ansatze for the classical background field have been proposed. The following 
table summarises what they will be referred to throughout the rest of the paper. 

the author therefore proposed a different ansatz, inspired by 't Hooft's multi- 
instanton form, that stays finite at the centre of the instantons, and dubbed it 
the ratio ansatz. It is given by 



In what follows we will refer to Re as the interactions or the ensemble generated 
by the ratio ansatz. We will compare the predictions from Re with those of the 
streamline ansatz S 48] and another 'hybrid' ratio-sum ansatz Rh 39]. This 
is summarised in Table [T] 

Phenomenological considerations have lead to the conclusion that the QCD 
vacuum consists of a dilute ensemble of instantons; a fact corroborated by lattice 
studies and self-consistency checks within the IILM. Diluteness and the localised 
nature of instantons render negligible contributions other than two-body inter- 



with O — 0\Oi- The formulas for an like-charged pairs follow trivially from the 
above. 

3.1. Gluonic interactions 

The complete classical gluonic interaction is given by the sum over all the 
possible pairings. It is clear from the structure of (|13l) that the colour degrees 
of freedom can be completely factorised out. After some lengthy algebra the 



4 A note on terminology: whenever we use the word interaction, we mean a quantity 'nor- 
malised' to the dilute gas, i.e. we subtract the dilute gas counterpart if the term naturally 
occurs in the exponential, as in the classical gluonic interactions, or we divide by the dilute 
gas counterpart if the interaction is a pre-exponential factor, as in the gluonic Jacobian or the 
Dirac determinant. 




1 + E* ^i{x,{yi,pi}) 



(12) 




(13) 



6 



15 



10 



R E : instanton-instanton 
R H : instanton-instanton 

R E : instanton-anti-instanton 

R H : instanton-anti-instanton 




3.5 



Figure 1: For instantons with equal sizes the interaction of Re agrees very well with Rh for 
oppositely charged instantons. There is slight discrepancy for like-charged instantons, in that 

the repulsion is a bit steeper in the Re case. (We have set p = \j p\ + p\ ■ ) 



result for the squared field strength can be written in the form 



pa pa 
P v fzis 



I + (TrO'O + (rjOtfi^^J + {f)Ori) PIJtpv I^ 



+ (riO t Ori) fipl , a J tip ^ cr + {■qOr]) apap (-qOri)p y p cr K ppy(J . (14) 

The different contributions are given in appendix | Appendix A| 

Factorising out the single instanton contributions and the coupling constant, 
the classical interaction between instantons is given by 



Sf 2 /S 



= V 12 = (S[A]/So - 2) 



S[A] = V 



pa pa 



(15) 
(16) 



where S[A] is the action of the background gauge fields and So = Sir/g 2 that of a 
single instanton. For equal sizes the agreement with Rh is very good, see Fig.[TJ 
However, for unequal sizes there are noticeable differences, see Fig. [2] The 
discrepancy follows from the functional dependence on the size parameters being 
of the form ' p\pi in Rh- As can be seen from the asymptotic behaviours, see 



appendix | Appendix A.2 the sizes enter rather in the combination \J p\ + pf, at 
least in the parameter regions of large and small separations; this is in agreement 
with [130. 



5 The authors of [Hjl have considered the sum ansatz; for large separations, however, every 
ansatz is equivalent to the sum ansatz. 
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Figure 2: For unequal size parameters, e.g. pi/p2 = 3 in this case, large differences start to 
become apparent. The reason is that the dependence on the sizes is more complicated than 
the functional form ^/p\P2 used in Ru- Note that the attractive well is much deeper in the 

Rh case which will eventually lead to a denser ensemble. (We have set p = , p\ + p|.) 



We will a adopt a couple of simplifications in the practical implementation 
that have been introduced in previous work. These are, on the one hand, the 
approximation of the high-frequency quantum interaction by an inverse running 
coupling constant evaluated at the scale of the mean instanton size and, on the 
other hand, the neglect of the Jacobian that introduces the collective coordinates 
and represents the low-frequency quantum interaction. 

In the single instanton case the high frequency quantum fluctuations lead 
to charge renormalisation and the coupling constant is replaced by the running 
coupling at the scale given by the instanton size [13], So(p) = &ir/g 2 (p). The 
same calculation has never been performed for a pair. In the original paper 
the interaction part of high-frequency quantum fluctuations have been estimated 
to be subdominant to the classical interactions . In that paper the authors argue 
then that the quantum interaction can be estimated by modulating the (total) 
classical interaction with the inverse running coupling constant, a slowly varying 
function of the background field, at the scale provided by the mean instanton 



size p. We will adopt the parametrisation put forth in [41], |39| that estimates 
the scale of the running coupling constant on a pair-by-pair basis and uses the 
geometrical mean of the sizes to set this scale. The full gluonic interaction is 
then given by 

S? 2 = S (^pTp7)V 12 ■ (17) 

The Jacobian interaction is positive by definition and can therefore be inter- 
preted as a repulsive (low-frequency) quantum interaction. A rough estimate of 
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the large distance behaviour suggests that the asymptotic power-law decay is 
0(l/i? 6 ), with R the separation between the pair. This is a faster decay than 
the well-known dipole-dipole interaction that follows from the classical action. 
For strong overlaps the Jacobian matrix will become approximately degenerate, 
and its determinant small, essentially because the matrix elements of the pair 
with the other instantons will be roughly equal. For complete degeneracy the 
singularity should be logarithmic because one singular value will tend to zero as 
the rank of the matrix decreases by one. The repulsion will thus be of the form 

In 7 sin s ro In ( R2 

with proportionality factor of order O(l) — O(10) because, as we argued, the 
degeneracy is due to one overlapping pair and should not get contributions 
from other instantons. In (|4.1j) we will discuss the small separation asymptotic 
behaviour for the ratio ansatz; the analytical expressions are given in appendix 
|Appcndix A.2.2| and we note that the singular behaviour is also repulsive and 
logarithmic, 

I*? ~ In (l + , (19) 

with a proportionality factor that is again of order O(l) — 0(10). In the in- 
termediate region it is harder to estimate the Jacobian interaction, but the 
logarithm should make its contribution subdominant. Also, the classical inter- 
action is boosted by the quantum contribution through charge renormalisation. 
We conclude that the Jacobian interaction is probably negligible compared to 
the classical interactions. 

Thus, the gluonic interactions we will use in this work will be given solely 
by the classical interaction. 

3.2. Quark Interactions 

The quark interaction arises from (jSJ, as is clear from our discussion in 
section [21 and is purely quantum mechanical. As for the gluonic interaction, 
some further approximations have been used in the literature; we will adopt 
these, albeit rephrased sightly differently. 

We will assume that the single instanton zero modes {£} form a functional 
orthonormal basis, i.e. we neglect contributions arising from non-vanishing over- 
laps among the £j. With this in mind, the finite dimensional low- frequency Dirac 
operator is then given by 

(p + m)ij = (d \p + m\Cj) = Pa + mSij . (20) 

To reiterate, we attribute the diagonal mass term to the requirement of or- 
thonormalitjffl rather than the degree of dilution of the instanton ensemble, e.g. 



6 Writing H = Hij\ipi) ® {*pj \ makes only sense if {ipi} forms an orthonormal system, given 
the scalar product (-|-). 
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371 ] . On the practical level this is irrelevant in as far as we recover the same 
determinantal interaction as used in previous works. 

The quark zero mode, in singular gauge and in the chiral representation, is 



given by [24 1 



f - 5sr^»ifE U, ) ■ (22) 

with <p Qa = e aa , normalised according to £12 = 1. Finally, Ui is the 3x3 
colour matrix describing the collective coordinates for the colour embedding; 
it is related to the adjoint representation by O ab = l/2Tr(t/T a £/V ), with 
U = U\U A - 

The Dirac operator, as defined above, is anti-hermitian. Eventually we need 
to diagonalise it, but, since readily available routines work with hermitian ma- 
trices, we display here the matrix elements of iJf). Within the ratio ansatz, the 
matrix elements Tja = J ^[ilnPn^A are as follows 

T IA = j^^^{Ur;)I,. (23) 

The concrete realisation of Ip is given in appendix | Appendix B| 

The rather large difference between Re and Rh, see top of Fig. [3l is due 
to the fact that the latter use a sum ansatz. The ratio ansatz was introduced 
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OR 



-25- 



/ 



R E : m u =0.1 m d =0.2 m s =0.4 
R H : m u =0.1 m d =0.2 m s =0.4 



/ 







2 



4 6 8 

Rip 



Figure 4: On the level of the effective interaction, the difference between Re and Rh is 
not as pronounced, i.e. the relative difference has decreased substantially. We can clearly 
see that light quark masses lead to a stronger attractive interaction between instantons and 
anti-instantons. Note that the relative difference between the ansatze Re and Rh does not 
seem to depend strongly on the quark masses. The instantons have been set up with equal 



to remove the unphysical divergence in the field strength; no such problem 
afflicts the overlap matrix elements. On top of that the quark determinant is 
a pre-exponential factor and as an effective interaction the extra logarithmic 
factor should make it rather insensitive to its exact form, see [H)]. Within 
our numerical framework, the full ratio ansatz does not produce any additional 
overhead and has the merit to be more consistent with the gluonic interactions. 
We have checked that upon neglect of the contributions special to the ratio 
ansatz, i.e. simplifying the overlaps so as to recover the sum ansatz, our results 
agree very well with those of Rh , apart form the aforementioned discrepancy in 
the instanton size parametrisation. As for the gluonic interactions, the colour 
matrices could again be completely factorised out. 

Note that the Dirac operator only connects instantons to anti-instantons 
due to the extra 7-matrix factor as compared to the mass operator, which van- 
ishes between instantons and anti-instantons. Therefore the quark fluctuation 
operator has the following form 



with T the Nj x Na matrix of overlaps Tja, and Nj {Na) is the number of 
instantons (anti-instantons); the O-matrices are Ni x Nj and Na x Na dimen- 
sional, respectively; finally, I is the identity operator on the quasi-zero mode 





(24) 
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space of dimension (Nj + Na) X (Nj + Na)- To diagonalise ip, it suffices to 
know the singular-value-decomposition of T. The left and right singular vectors, 
vp L and ip R , are defined by 



(25) 
(26) 



The singular eigenvalues A„ are always positive. The kernel of the Dirac operator 
is spanned by the A = singular eigenvectors, ip K , of either T ox T\ depending 
on whether N] < Na or Nj > Na- We can then construct the eigenvalue 
decomposition of the Dirac operator. The non-zero eigenvalue part has the 
following eigensystem 



— A„, 



€{1, 



Finally, the kernel is spanned by the eigensystem 
0. 







0. 







ne{l,...,JV A -JVj} 
nG{l,...,JV>-^} 



a(Arj,7V A )} 

N T <N A , 
N!> N A - 



(27) 



(28) 



Note that the non-zero eigenvalues come in pairs. Together with the zero 
eigenvalues, the determinant of the Dirac operator can be written as 



mm(iVi,AU) 

det(ip) = m IQI J| (m 2 + A„) . 



(29) 



with Q = Nj — Na the topological charge. If we are only interested in the deter- 
minant, and not so much in the eigensystem, this can be put in the equivalent 
form 

{Ql j det(TTt+™ 2 ), Q<0 

det(Ttr + m 2 ), Q > ' 



(30) 



Upon placing this term into the exponential, the normalised determinant of 
quark zero mode overlaps leads to an effective interaction. The normalisation 
consists of dividing (I3U1) by m I+ A . The quark interaction is thus given by 



q lndet(TTt + to 2 ) -iV/lnm 2 , Q<0 

N f " 2^ i In det(TtT + ml) - N A In m 2 , Q > 



(31) 



with Nf the number of active quark flavours. Note that the quark interaction 
is always attractive. This follows from the fact that we can write the overlap 
matrix for each flavour as I + and this form makes it explicit that the 
determinant is bounded from below by unity because the smallest eigenvalue is 
easily seen to satisfy A m ; n > 1. 



12 



This exhausts the interactions in the IILM because the fluctuation operator 
of the ghost part is positive definite and its lack of zero modes prevents the 
construction of the low frequency part of the spectrum within the moduli-space 
approximation. We are thus left with the high-frequency part which, as in the 
other cases, is assumed to factorise and cannot lead to interactions. 

4. Numerical Implementation 

4-1- Interpolation and asymptotic matching 

The decoupling of the colour degrees of freedom is a computational benefit: 
by using global 50(4) transformations, without loss of generality, we place the 
first instanton at the origin and the second along the z-direction. The initial 
orientational dependence is then factored out of the integrand and combines 
with the colour matrices as in [l3j . These integrations are too time consuming 
to perform during actual simulations; instead, they are computed beforehand 
to fill interpolation tables that are, in turn, used during the simulations. The 
interpolation grid is three-dimensional, and depends on p%, p2 and R — \x\ — X2\- 
For numerical stability we choose to use simple linear interpolation. 

A uniform grid can, of course, only extend over a finite region and we must 
decide which portion of the parameter space to cover. We took the single in- 
stanton moduli-space measure as a guide for the size grid because, suitably 
normalised, it can be interpreted as a probability density. We choose the lower 
limit, p m j n « JjjA, to be a fairly small quantilqj- Here, A is the scale at which 
QCD starts to become strongly coupled. The upper limit is set to p max = A. 
Larger instantons cannot be treated consistently in the IILM because it uses 
perturbation theory, which breaks down below A. 

We believe that these choices cover the relevant parameter space, and we 
sample the sizes from the interval [praxm Pmax]- As a consistency check we moni- 
tored the actual size distribution and did not find any evidence for a significant 
weight at the edges of the sample interval. We therefore conclude that this 
procedure is well-defined. 

The classical gluonic interaction in the ratio ansatz suffers from gauge sin- 
gularities that prevent us from extending the grid down to vanishingly small 
instanton separations, R —¥ 0. The opposite limit, R — > oo, cannot be covered 
either unless we use a non- uniform measure on R + . In principle this would seem 
like the most elegant approach, however, it is not feasible practically because the 
numerical integration becomes inaccurate at larger separations; the only rem- 
edy would be to set very small error tolerances for the numerical integrations, 
but that is computationally prohibitive. Therefore, we decided to use matching 
formulas for both the large and small separation regimes. 

The rationale is not to derive accurate formulas in absolute terms but to get 
the absolute value from the interpolation results at a matching point R rn . The 



7 It corresponds to less than the millionth quantile. 



13 



Figure 5: The instantons Ii and I2 are so far apart that, within the shaded region that give 
the dominant contribution to the field strength of each, the partner's field strength is roughly 
constant and fixed at x M — i? M — We can then safely extend the integration region 
to be all of M 4 , with a negligible error due to the rather strong localisation of the individual 
instantons. 



matching formulas are thus to be understood as accurate in a relative sense, i.e. 
the asymptotic interactions / asy should behave, asymptotically, like the exact 
numerical interactions / ox . This ensures that we reproduce the correct fall-off 
or singularity behaviour. Thus, we compute the interactions according to 



/ex(-^m) 
/asy \Rm ) 



f(R) = UiR) i cxy Z m \ , (32) 



whenever they fall out of the grid. Since the localisation of the instantons is set 
by the sizes, it is natural for the matching point to be proportional to the for- 
mer. Eventually, the exact proportionality factor follows from an 'optimisation' 
procedure, given that we aim for the interpolated interactions to be correct at 
the one percent level. 

The full gluonic interaction consists of different pieces that are added to- 
gether, (fl4|) . We could use (|32|) for these subinter actions term by term but 
it turns out that such a matching is numerically rather unstable. Thus, even 
though we are only interested in asymptotic relations, we need a systematic 
procedure that insures that the different asymptotic subintcractions are added 
up with the correct magnitude relative to each other. 

For the large separation case we want the instantons to be so far apart from 
each other that within the region in which the field strength for I\ is stron, 
the field strength of I2 hardly changes: we can approximate — w — R 
Since the field strength is negligible at and beyond R^, we can safely extend 
the integration region to cover all of M 4 . The field strength of I2 behaves as a 
constant, and we can use the rather simple rational expression for the interaction 
in terms of the 't Hooft potential to find exact results. We add to this the 
analogous contribution from I± Jg. The configuration is illustrated in Fig.[SJ 

We shall call this the zeroth order approximation, and it is clear that, to this 
order, terms odd in derivatives of Ii will vanish due to 0(4) symmetry. However, 



We use a translation to place I\ at the origin. 
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Proper small separation asymptotic integration down to R (I) 




Large separation asymptotic integration up to ^ (II) 



Figure 6: The instantons I\ and I2 are strongly overlapping. We approximate the integral 
by, first, integrating over I\ keeping I2 fixed at _R M , as in the large separation case but with 
upper limit R/2; to this we add the analogous contribution from 1%. Secondly, the, possibly, 
singular behaviour is picked up by integrating from infinity down to R, and approximating 
the arguments to be x^ — Rfi/2 x M and x M + R/j/2 w x^ respectively. 



it turns out numerically that, upon combining all the different terms from (I14[) , 
the zeroth order terms are sufficient. In particular, no non-integrable terms 



are present; as we will demonstrate in [54lp. the finite temperature interactions 



do not include terms such as ([T]) that prohibit the thermodynamic limit. Our 
formulas are given in appendix |Appendix A.2.1] 

In principle, we can compute the neglected terms by going to first order, i.e. 
g(x — R) m g{—R) + x^d^g(-R), or beyond. Such higher order contributions 
will typically no longer converge on M 4 . It seems natural to cut them off at 
i?, and this will generally lead to logarithms, ln(l + R 2 /{p\ + p$)), together 
with rational functions. However, in contrast to the fitting formulas of (39l |. 
the Taylor expansion of our asymptotic formulas produce only power-law like 
decays for large separations; in addition they fall off more strongly than the 
zeroth order terms and thus will not produce non-integrable interactions either. 
We have thus achieved our goal of deriving interactions that allow us to study 
the thermodynamic limit of the IILM. 

We now turn to the case of asymptotically small separations. A typical 
situation is depicted in Fig. [5] The rationale is to split the integration into 2 
regions. 



9 At finite temperature, the present framework remains virtually unchanged. The only 
modifications are that the 't Hooft potential, II, changes and that the integration region 
becomes S 1 X I 3 . 



15 



10 
10" 1 
10" 1 
10" 1 



10 v 



Gluonic interaction (approximation) 
Gluonic interaction (exact) 
Quark interaction (approximation) 
Quark interaction (exact) 



10- 



10- 1 



10 



Figure 7: The gluonic interaction is well approximated by the combination of interpolation 
and asymptotic matching. The quark overlap is very poorly approximated by the zeroth order 
small separation asymptotic formula; it tends to zero with too high a power as compared to the 
exact result. The correct behaviour can in principle be obtained from higher orders, and we've 
estimated that the second order contribution will suffice. In practice, the quark interaction in 
this region is completely irrelevant as compared to the gluonic interaction. 



I The far-field region beyond both centres, placed symmetrically around the 

origin; we approximate the arguments by ± R^/2 sa ±x M . 
II The region around each instanton up to R/2, with R the pair separation. 
We integrate around ±R^/2 keeping the arguments of the partner instan- 
ton fixed at =F R^/2 « =Fi? M /2. This is similar to the large separation 
case, but here we only integrate up to R/2. 

Region I accounts for possible singularities. After adding up all the different 
subintcractions, the singularities from region I dominate the total interaction. 
Since we need the region II approximations anyway in the large separation 
case, it does not represent any extra overhead to use them as well in the small 
separation limit. 

In Figs. [7] and [5] we plot the exact and approximate result for the gluonic, 
quark and total interaction. Note that some subinteractions in the gluonic sector 
are poorly approximated by the zeroth order asymptotic matching formulas in 
region I. However, those terms that do exhibit singularities completely dominate, 
and the total gluonic interaction is well approximated for all separations. The 
quark overlap consists of just one term, which is not well described in region I. 
We estimated that the correct asymptotic behaviour can be obtained at second 
order. However, this is unnecessary since the quark interaction is bounded 
in that region and the gluonic interaction completely dominates. Combining 
the gluonic and fermionic interactions, we see that the total pair interaction is 



1G 




Figure 8: The total pair interaction is accurate on the one percent level. Note that the error 
in region I from the quark interaction is completely negligible. The spike in the bottom plot 
is due to the sensitivity to zero crossings. 
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accurate on the one percent level. 

Note that we use interpolation up to very strong overlaps; also, instantons 
rarely enter region I during simulations because it makes up a very small part 
of the total volume box. In any case, region I overlaps will almost certainly 
be rejected in Monte Carlo moves, and therefore quantities that are computed 
solely from the quark interaction, as for instance the quark condensate, are very 
insensitive to large errors in region I. Ultimately, this is the reason why we chose 
the interpolation grid to cover such strong overlaps. 

4.2. Monte Carlo 

Previous studies, lattice results and phenomenology indicate that the instan- 
ton ensemble is fairly dilute. Therefore, we organise the partition function into 
a dilute gas measure times the exponential of interactions, 



Z = 



E 

Ni,Na 

OO 

E 



1 1 

ivji Na\ 



1 1 



Ni,N A 



S" = J2 S ; 



T W7\ Ni ' Naj 



TV/! N A 



(33) 

(34) 
(35) 



where Sf^ is given in (|17l) and S q is given in (|31l) . We follow [3^] and use one- 
loop accuracy for the charge renormalisation factor that modulates the classical 
gluonic interaction, i.e. So (-v/plpi) — ► P1P2), with j3\ given below in (|4"Tj) . 
Although not really consistent, the single instanton density is given at two-loop 
in order to replace the pre-exponential bare by running coupling constants [l3j ; 



the former were induced by the transformation to collective coordinates, 
two-loop single instanton measure is then given by 



d{p)=dl{p)dl{p) N s , 
d 9 (p)=C Nc p~ 5 f3 1 (p) 2N ^exp 
1 



-ft(p) 



2N r 



b'\ b' ln/3i{p) 



2b J 2b 0i(p) 



c/q(p) = mpexp 



• In mp 



3 In mp + 2a — (6a + 2/3)(mp) 2 + 2A 1 (mp) A - 2A 2 (mp) 6 * 
1 - 3(mp) 2 + B x {mpY + B 2 {mpf + B 3 (mp) 8 , 



The 

(36) 
(37) 

(38) 
(39) 
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For the quark term, dp, we use the generalisation of 't Hooft's j47| result valid 
for arbitrary mass [11| . The different terms in d 9 are given by 

0.466 e" 1679 ^ fAn . 

Cn < ~ { Nc -mN c -2)\ ' (40) 

ft (p) = -b ln(pA) , b = y N c - 2 -N s , (41) 

Note that the above has been derived in Pauli-Villars regularisation. 

Being an interacting many-body system, the partition function cannot be 
evaluated analytically, and we choose Monte Carlo methods to cope with it 
numerically. More precisely, we will use the Metropolis algorithm to sample 
the important integration regions of the partition function. This is, of course, 
all well known, but it seems appropriate to introduce the, possibly less known, 
Monte Carlo moves corresponding to insertion and deletion of instantons needed 
for grand canonical simulations. 

Following the usual strategy of imposing detailed balance, the simplest in- 
sertion/deletion algorithm consists of randomly placing an instanton in the box 
and randomly selecting an instanton to be removed. Imposing detailed balance 
and considering the case of an instanton, we arrive at 

yPN t ,N A -ANi,Nj+l = N + 1 P < Ni+1,N a -Ani+1,Ni ■ (43) 

As usual, p ^ N — Zn i n a I^\ is the probability to be in the state {Ni,Na}- 
The acceptance probability Aij is implicitly defined through (|43p. and the 
Metropolis algorithm defines it to have the following form, [19j |. 

•AiVj.i^+i = min(l,^l), (44) 
A Ni +i,n i = min(l,^ 1 ). (45) 
Plugging this into (|43| we finally arrive at 

-4= / + /i +1 ^ - (46) 
Ni + l Pnj.n a 

The difference to ordinary Monte Carlo moves, as used in the canonical 
ensembld"! is that the proposal probabilities do not cancel and the transition 



"Note that we neglect the factorial terms in the definition of the equilibrium probability 
density p eq because they are an artifact as far as the measure is concerned. They have been 
introduced to render the integration volume simple, i.e. the product of the single instanton 
moduli-spaces M 1 . During the integration process all the permutations of a given set of 
coordinates are generated, but, since the instantons are indistinguishable, they really corre- 
spond to only one configuration. To correct for this overcounting, we then have to divide by 
a factor of Njl. The important point is that for the transition probabilities these factorial 
factors are irrelevant. 

J1 That is, updates for the positions, sizes and colour orientations 
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5c 


See 


(N) 




(Q 2 ) 


& 


(Shit) 


£int 


0.5 


0.5 


101.4(6) 


2500 


2.4(1) 


200 


-5.005(5) 


1200 


0.6 


0.4 


100.8(6) 


2100 


2.30(5) 


130 


-5.0010(5) 


1400 


0.8 


0.2 


102.6(6) 


4000 


2.35(4) 


130 


-5.010(2) 


2400 


0.9 


0.1 


102.1(6) 


5000 


2.5(1) 


270 


-5.016(3) 


3700 



Table 2: The sample size is roughly equivalent for each set, with 200 independent configura- 
tions generated according to the autocorrelation time £jy Considering some bulk properties, 
we see that the sampling does not really depend on the a-priori-probabilities Si . Even though 
the autocorrelation times are only rough estimates, we will take the data at face value and 
choose 8c = 0.6 and Sqc = 0.4 for the remaining simulations. 

matrix is not symmetric. In this specific case, the proposal probability for an 
insertion is V l p s = 1/V, corresponding to the probability to place the instanton 
randomly within the box, whereas the proposal probability for a deletion is 
Vp cl = l/(Nj + 1), corresponding to the probability to select an instanton 
among the Nj + 1 available. 

When we perform the standard updates, it is easy to monitor the acceptance 
rates and tune the the proposal probabilities to achieve good rates, i.e. 50% say. 
For the move described by (|46| we do not have a parameter to tune though. At 
T = 0, this is not really a big issue because it turns out that the acceptance 
rate is ss 0.4, even for the rather small quark masses that we will use. This is 
still acceptable and does not really justify the overhead of more sophisticated 
update algorithms. 

Also, note that such an insertion/deletion step is a fairly large change as 
compared to the normal coordinate updates, and so these grand canonical moves 
actually help to sweep through phase space more quickly. 

Finally, we need to decide how many grand canonical moves we perform per 
coordinate update, that is we need to fix the a-priori-probabilities 5c and 5gc- 
We found that, for T = 0, the ensemble is not sensitive at all to this parameter, 
see Table [2] Since we will ultimately be interested in computing the topological 
susceptibility, we will aim to achieve low autocorrelation times for the instanton 
number, N = Nj + Na, and topological charge, Q — Ni — Na, i.e. we will 
perform rather more insertion/deletion moves than less. In practice we perform 
canonical moves only 60% of the time. 

4-3. Fermionic determinant 

As mentioned in the introduction, we want to study the IILM for 'physical' 
quark masses. In that case, we must make sure that the simulation box is large 
to be insensitive to finite size effects. In the lattice community it is common 
practice to use a box length that corresponds to 4— 5 times the wavelength of the 
lightest propagating degree of freedom, which is the pion. In practice, we want to 
circumvent the need for extremely large boxes by studying the thermodynamic 
limit, V — > oo. 
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As compared to fitting formulas, our combination of interpolation and asymp- 
totic matching results in a rather substantial computational overhead. This is 
particularly so in the quenched case. For unquenched simulations the situation 
is less drastic as the computationally most demanding part is the evaluation of 
the determinant and/or the determination of the eigensystem of the Dirac op- 
erator. Increasing the simulation box, i.e. increasing the number of instantons, 
this becomes the bottleneck to large volume simulations. 

The Monte Carlo changes are, however, of a rather simple form, changing 
only one column of the overlap matrix T at a time. We can therefore use de- 
composition update techniques to reduce the complexity from 0(N 3 ) to 0(N 2 ). 

For the update step we only need to evaluate the determinant (j3"0")) . Given 
the fact that m 2 + TT', respectively m 2 + T^T, is a positive definite hermitian 
matrix, the fastest evaluation will be achieve by using the Cholesky decompo- 
sition. An added bonus is that the Cholesky decomposition and its algorithm 
are known to be very stable. 

Focusing on M 2 = m 2 + TT^ = LDL\ an update V = T + AT can be 
written as two rank 1 updates for M 2 , of the form 

A/' 2 = M 2 + $$ f - ^ , (47) 



with $, \fr vectors. Details are given in appendix Appendix C where we also 
discuss more efficient ways to deal with adding and removing instantons, and the 
corresponding updates. The Cholesky decomposition can be updated efficiently 
when it only changes by rank 1 matrices, that is transformations of the form 

M' 2 = L'D'L ^ = M 2 + azz^ = L(D + aww^tf , (48) 

where Lw = z. The algorithms then compute the decomposition of D + aww* = 
LDL\ which can be achieved in 0(N 2 ) because D is diagonal. Furthermore, 
the matrix L has a special form which allows an efficient matrix multiplication, 
V = LL, in 0(N 2 ). Details can be found in [22j. The algorithm we use in 
practice is known to be unstable for downgrading, a < 0, unless the resulting 
matrix, M' 2 , is known to be positive definite. Since upgrading, a > 0, is always 
stable, it is important to perform the two consecutive updates in the order given 
by (E]). 

In general we will be performing grand canonical simulations, and need to 
keep track of two decompositions, one for m 2 + TTt and one for to 2 + T^T. 
Furthermore, we deal with 3 active quarks so that each Monte Carlo update 
entails 2 ■ 2 • 3 = 12 rank 1 updates. We find that for an ensemble with 100 
instantons and 100 anti-instantons we achieve a computational gain of a factor 
of 2 as compared to the full Cholesky decomposition. 

5. Different Ensembles 

To be predictive, the IILM should not depend too sensitively on the chosen 
ansatz. Given the fact that, for instance, the streamline and the ratio ansatz 
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have quite different functional forms for overlapping pairs, the insensitivity of 
the model to specific background ansatze can only be determined a posteriori. 
On a heuristic, level we expect insensitivity to emerge if the ensemble stabilises 
in a rather dilute form so that the precise functional form of the repulsion 
is irrelevant. The large separation limit is a priori unproblematic because all 
ansatze are constructed such that, asymptotically, they approach the simple 
sum ansatz, i.e. A = A\ + A%, with A the gauge field. 

First, we will frame our discussion on the pair interactions. The total effec- 
tive interaction of a pair of oppositely charged partners is given by 

S12 = S (^T 2 )V l2 - ^ In ( \ T ^\ 2 + m l \ . (49) 

n=l \ m n / 

Identically charged pairs only feel the gluonic interaction as Tjj = Taa — 
0. As expected, the ratio and streamline ansatz are markedly different only 
for strongly overlapping pairs, see Fig. [SJ where strongly overlapping pairs are 
characterised by R < yj ' p\ + p\- 

In the quenched case, we notice that the Re ansatz has a higher absolute 
minimum as compared to the Rh ansatz, occurring roughly for the same sepa- 
ration. So we expect the ensemble to become slightly more dilute because it will 
not be as favourable, energetically, for instantons to come close. For unequal 
sizes, however, the repulsion is weaker in the Re case which would favour a 
denser ensemble, as less volume is excluded. The streamline ansatz will lead 
to a substantially more dilute system because the core repulsion is broader, 
excluding more volume for the instantons to move through. 

In the unquenched case, the difference in the absolute interaction strength is 
much more pronounced between Re and Rh- We therefore expect that the Re 
ansatz should be quite a bit more dilute as compared to the Rh ansatz. Con- 
sidering that the streamline ansatz has a deeper minimum than the Re ansatz, 
the former will favour instantons to come closer. However, it has more excluded 
volume. Both trends work in opposite directions, and there is a possibility that 
they lead to roughly identical ensembles, at least on the level of the instanton 
density. 

We will address these issues in more detail by perfor ming canonical simula- 
tions and minimising the free energy. This follows closely [39J , and also serves to 
validate our code against their results. Note that the simulations are performed 
in the topologically trivial sector, for which Nj — Na, i-e. the topological charge 
Q = 0. In Fig. [10] we plot the free energy F = — \nZ/V against the instanton 
density n = (Nj + Na)/V. As expected from our considerations of the pair 
interactions, in the quenched case the Re and Rh ansatze are only slightly 
different, with Re lea ding to a slightly denser ensemble. Also, the interactions 
stored in that ensemblqlj are a bit lower, again as could be anticipated from the 



1 The difference in the free energies is directly related to the difference of the interaction 
per instanton 
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Figure 9: Top: Most attractive colour orientation. Bottom: Random colour orientation. For 
both graphs the instantons have different size parameters. We expect the Rh ensemble to 
be denser than the Re ensemble because the attraction well is deeper, whereas the excluded 
volume due to repulsion is not that different. Along the same lines, the S ansatz should lead 
to a rather more dilute system in the quenched case. For unquenched simulations, the deeper 
attraction well, and the steeper and broader repulsion, of the S interactions might lead to an 
ensemble roughly equivalent to the Re one. 
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m u 


md 


m s 


Mi 


0.1 


0.1 


0.7 


M 2 


0.05 


0.05 


0.3 


M 3 


0.012 


0.022 


0.44 



Table 3: We use three different sets of quark masses and investigate how the instanton liquid 
depends on them. 

pair interactions. The S ansatz leads to a much more dilute instanton ensemble, 
and our data reproduces well that of 39] . 

Ultimately, we will be interested in smaller quark masses. It is clear from 
(j4"5)) and Fig. 2] that smaller masses increase the quark interaction strength as 
compared to the gluonic counterpart, which stays constant and is responsible for 
the core repulsion. From a purely energetic point of view, smaller quark masses 
should then lead to denser ensembles. However, we clearly see in Fig. [10] that 
the ensembles become more dilute. The reason is that the small quark masses 
enter the instanton size distribution; in turn, the density, in the dilute gas limit, 
is entirely set by the size distribution, i.e. n — 2 J dpd(p). Bringing it into the 
action, we can interpret the size distribution as the energy cost needed to insert 
an instanton into the box. This is a well-known fact, namely that small quark 
masses suppress instanton contributions to the QCD vacuum because the differ- 
ent topological vacua become equivalent in the limit of vanishing quark masses; 
phrased differently, the energy barrier has disappeared, and only field configu- 
rations with topological charge Q = survive. In the dilute gas approximation 
this leads to the disappearance of instantons altogether. As Fig. 1101 shows, this 
is not true for an interacting instanton ensemble, where the instanton density 
converges to a finite limit as the quark mass is lowered^. The results from 
Fig. [10] also show that the Re ansatz generates an ensemble that differs more 
and more from the Rh ansatz, as was anticipated from our considerations of 
the pair interactions. We can also clearly see that the Re ensemble does not 
converge to the S ensemble. 

So far we have framed the discussion essentially in terms of the instanton 
density. To investigate the similarities and differences in more detail, we will 
look at a few bulk properties and their dependence on the different ansatze in 
the thermodynamic limit and the grand canonical ensemble. We clearly see 
how the density decreases with the quark masses, but approaches a finite limit, 
Fig.QT] As we have discussed before, the quark masses will suppress fluctuations 
to inequivalent topological sectors and in the limit of vanishing masses only the 
trivial Q — sector will survive. The fluctuations between topological sectors 
are encoded in the topological susceptibility \ — (Q 2 ) /V, which vanishes with 
the quark masses, see Fig. Q2J Both the instanton number and the topological 



13 Remember that the simulations take place in the topologically trivial sector, i.e. Nj = Nj\ 
or Q = 0. 
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Figure 10: The different quark masses are summarised in Table \3\ The simulations were 
performed in the topologically trivial sector, i.e. Ni = Na- For the quenched and the Mi 
case we fixed N = 64 as in |3Sl | - The other two unquenched simulations have N = 200. We 
clearly see that small quark masses suppress instanton contributions to the QCD vacuum, 
but also that there exists a finite limit for the instanton density as the quark masses vanish; 
this is in contrast to the dilute gas approximation which suppresses instanton contributions 
completely for zero quark masses. For unquenched simulations the free energy for the Re 
ensemble roughly agrees with that of the S ensemble, although the equilibrium densities are 
rather different. Still, the approximate equality between the free energies might be interpreted 
as evidence of an approximate equivalence between both ensembles for bulk properties, e.g. 
equivalent pressure, since it is directly related to the free energy. However, the Re liquid does 
not seem to converge towards the S ensemble as we lower the quark masses. 
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Figure 11: As anticipated from the canonical study, the instanton number for both ratio 
ansatze is very similar in the quenched case. There seems to exist a finite limit for the 
instanton density as the quark masses vanish, and instantons will be present in the QCD 
vacuum even in the chiral limit; this is in sharp contrast to dilute gas approximations. 

charge fluctuations exhibit a nice scaling with the volume. 

We will now turn to an intensive quantities, the mean instanton size p Fig. 1131 
For the quenched and the first unquenched, Mi, simulations, which have been 
calibrated to achieve N ~ 64 for the smallest volume as in the canonical simu- 
lations, we see that the simulation boxes are not large enough, even though the 
density and the charge fluctuations seem to suggest otherwise, i.e. display good 
scaling with V. For the other two unquenched simulations, tuned to N m 200 
for the smallest volume, we have reached volume sizes large enough to perform 
a thermodynamic limit. It is worth noticing that the mean instanton size is a 
rather robust quantity, and does neither depend strongly on the ansatz nor on 
the quark masses. This makes it a good quantity to use when comparing data 
from different ensembles, e.g. Q where the authors establish that the scale at 
which the IILM is operating is given by the inverse of the mean instanton size. 
Another intensive quantity, the interaction per instanton, is less sensitive to 
finite size effects, see Fig. QjJ The data shows that the weaker repulsion of the 
Re ansatz as compared to the Rr ansatz dominates over the deeper attractive 
well of the latter. Therefore, the total interaction in the Re ensemble is slightly 
lower, leading to a denser system. The stronger repulsion for the S interactions 
leads to more excluded volume; this, in turn, leads to lower interactions and 
a more dilute ensemble. These conclusions are in agreement with the direct 
measurement of the instanton density Fig. 111! 

We are mostly interested in unquenched results, and the following comments 
relate this sector. From the data of Table HI we can infer that results for x 1 ^ 4 
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Figure 12: The topological susceptibility, the slope of the graphs, is very sensitive to the 
quark masses. It is screened by small quark masses and will vanish in the chiral limit. This is 
expected as QCD with massless quarks does not have topologically inequivalent vacua; in this 
case the so-called 6 parameter is not physical and can be rotated away by a chiral rotation of 
the quark fields. See also [4d | for another work on the topological susceptibility in the IILM. 
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Figure 13: For the quenched and the Mi simulations, which have been fixed to N fa 64 for the 
smallest volume, the simulation boxes are still a too small, as can be seen by the systematic 
drift. For the other two unquenched simulations (N 200 for the smallest volume) we are 
much closer to the thermodynamic limit, although there are still systematic deviations for 
the Rh ansatz. In any case, the different ansatze give rather similar results. Also the mean 
instanton size approaches a unique limit for small quark masses. 
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Figure 14: As anticipated from the pair interaction considerations, the interactions are very 
similar in the quenched sector for Re and Rh- For full simulations, the differences between 
the ansatze stay constant as the quark masses vary, with Rjj leading to the strongest attractive 
interactions, as was expected. Note that the pair interactions are less sensitive to finite size 
effects compared to the mean instanton size. 





n 


X 


(P) 


(S g + S<*)/V 


Quenched 


0.567(1) [Re] 
0.532(1) [R H ] 
0.282(1) [S] 


0.46(3) 
0.86(5) 
0.24(1) 


0.6837(2) 
0.6850(1) 
0.5631(1) 


1.420(1) 
1.455(1) 
1.122(1) 


Mi 


0.288(2) 
0.370(2) 
0.163(1) 


0.041(1) 
0.0369(8) 
0.0196(6) 


0.6662(2) 
0.6581(3) 
0.6331(3) 


-1.648(1) 
-2.311(2) 
-1.818(3) 


M 2 


0.1660(7) 
0.259(1) 
0.1255(8) 


0.0136(4) 
0.0126(2) 
0.0075(2) 


0.6757(2) 
0.6615(2) 
0.6510(2) 


-2.996(1) 
-3.945(2) 
-3.472(2) 


M 3 


0.1686(8) 
0.265(1) 
0.1269(5) 


0.00440(7) 
0.00403(7) 
0.00230(5) 


0.6744(2) 
0.6606(2) 
0.6511(2) 


-4.954(1) 
-5.902(2) 
-5.407(2) 



Table 4: Thermodynamic extrapolations for the instanton density, the topological suscepti- 
bility, the mean instanton size and the mean interaction. The data has been obtained from 
Figs . HT1H21H51 and [14] respectively. 
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have a 12% systematic ansatz dependence. The topological susceptibility is 
surprisingly similar for the Re and Rr ansatz in the unquenched sector. The 
mean instanton size is indeed a rather robust quantity and only affected on the 
3% level by these systematics. Also, note that the mass dependence on (p) is 
rather small, with differences not larger than 5%. The instanton interactions 
and the n 1 / 4 agree within 20%, and the latter converges to a fixed limit as the 
quark masses vanish. 



6. Fixing parameters 

6.1. Quenched case 

In the quenched case, the IILM has only one freely adjustable constant, the 
lambda parameter A. We need one observable, from the lattice say, to fix it. 
Different approaches can be chosen. In the early works, A was determined by 
fixing the instanton density to 1 fm~ 4 at T = 0. To compare this with the lattice 
is not straightforward as the classical instanton content is convoluted with the 
quantum mechanical fluctuations. With the discovery of the KvBLL calorons, 
there is a renewed interest in studying the topological structures on the lattice, 
see for instance Since the topological susceptibility is well measured on the 
lattice and is easily accessible within the IILM, it is a natural candidate. The 
lambda parameter is then given by 



A=4 /Jfl«*_ - ( 5 o) 

V XIILM 

We will use X\it = 193 MeV, [12]. The topological susceptibility in the IILM is 
extracted from Fig.[T5]by using the definition 

Xtap=limM (51) 

V — ►oo V 

This yields A = 234(1) MeV. The error is purely statistical. The instanton 
density turns out to be n — 0.543 A 4 = 1.02(2) fm -4 , fairly close to the usually 
quoted phenomenological value of n = 1 fm~ 4 . We find that even for these 
larger volumes the mean instanton size is still evolving towards lower values, 
as in Fig. [T5J The largest volume then leads to the upper bound p < 0.57 fm. 
Using a simple fit to p = p^ + aV~ ' 25 to extrapolate to the asymptotic value, 
we find poo ~ 0.53 fm; this is rather large compared to the phenomenological 
value of p ~ 0.33 fm. 

To estimate the systematic error due to the dependence on the ansatz, we 
will use the data from Table |U The fact that we take a fourth root reduces the 
rather large differences in xiilm to about 15% for A, i.e. A = 234(35) MeV. Our 
value has been obtain through simulations in PV regularisation. To compare 
it with lattice data, we will convert it to the MS scheme, [26|, A^g/Apy = 
exp(— 1/22). This gives A^g- = 224(33) and compares well with the lattice 
result A M g = 259(20) j23|. 
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Figure 15: The fluctuations of the topological charge (Q 2 ) show a nice linear dependence with 
the volume box V, as it should be for an extensive quantity. From this we infer the topological 

susceptibility Xtop = limy-nx, ^jr 1 - ■ 



6.2. Unquenched case 

We want to use realistic quark masses. These are fairly small, and one must 
worry whether such light degrees of freedom will fit into the simulation box. 
The usual approach, used in the lattice community and also in work on the 
IILM is to compute the pion mass from a set of unphysical quarks and to 
fix the volume box such that Lm^ > 5; chiral perturbation theory can then be 
used to extrapolate to physical masses. Ultimately, the lattice wants to test 
the predictions of chiral perturbation theory as well, and in recent years, the 
computing power and, most importantly, the algorithms have improved to such 
an extent that physical quark mass simulations are becoming feasible; however, 
these are still immensely costly simulations, and 2 + 1 flavour simulations were 
rare until recently. 

We follow a rather more modest rationale by simply demanding that the 
quark mass be at least so small as to be comparable to the lowest eigenvalue 
of the Dirac operator, (A m ; n ), see Fig. [17] This sets the smallest box we use 
in our simulations. We then use ever larger volumes and extrapolate to the 
thermodynamic limit. 

In Q the lambda parametei0 is fixed by computing the meson and nucleon 
masses, through current correlators of the interpolating fields and their asymp- 



14 Actually, the authors fixed the mean instanton size. But it is trivial to relate the latter 
to the lambda parameter. 
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totic spatial decay, and by comparing them with the available lattice data. This 
study established that the IILM is compatible with the predictions of chiral 
perturbation theory. We will take this for granted in what follows. 

In order to fix A, we could still use the topological susceptibility as it is rou- 
tinely measured on the lattice. However, the topological susceptibility depends 
strongly on the quark masses, see Fig. Q21 We can get rid of the mass depen- 
dence by using chiral perturbation theory and computing the chiral condensate 
(990- The chiral condensate has been studied within chiral perturbation the- 
ory and, more recently, it has been precisely determined on the lattice @, H, HI- 
We will take it to be (qq)^^ = 2 GeV) = 250 MeV. 

To extract the chiral condensate from the IILM, we will use the procedure 
adopted in 0, 0, @ : we compute the topological susceptibility for different sets 
of quark masses and extrapolate to the chiral limit, m, — > 0. The condensate 
can then be determined by chiral perturbation theory [33] , 

X = m cS (qq) + O(m 2 ), (52) 



The chiral condensate has an anomalous dimension and, therefore, depends 
on the scale. Furthermore, the IILM is set up with a PV regulator, whereas 
the quoted result is computed in dimensional regularisation. It is well known 
that within an unphysical renormalisation scheme such as MS^l the results de- 
pend on the regulator (for unphysical quantities like masses, coupling constants 
and amplitudes). We therefore need to compute the finite counterterms that 
relate the PV to the MS regularised results. Deferring the details to appendix 
|Appendix D] we find that (qq)„ v (fi = 2 GeV) 244 MeV. This is a one-loop 
result. The two-loop correction can be estimated very roughly to be on the 10% 
level as is typical for computations around the scale of \i = 2 Ge\0. 

We will define the scale of the IILM by /i\ = A/p, as suggested in Q, and 
determine A from the self-consistency equation 

(q q f V M=A 3 (qq)^ M , (53) 

where we run the chiral condensate (qq)Q V (/i) at one loop. To that order, there 
is no difference between schemes and we can use the MS results, e.g. [49j | . 




5 To reiterate, we implicitly rely on the fact that the IILM is describing well the chiral 
properties of QCD, as has been checked in numerous studies, the most convincing being (s|. 

16 't Hooft's computation of the one- loop instanton measure, using Pauli-Villars regular- 
isation, is also unphysical because, instead of poles, logarithms of the regulator mass are 
subtracted. 

17 Strictly speaking, we should use the two-loop result because the simulations in the IILM 
have been obtained using the two-loop improved instanton measure. However, Pauli-Villars 
regularisation is not straightforward for non-Abelian gauge theories beyond the one-loop level, 
and we do not have the expertise to embark on this endeavour. In any case, the difference 
should still be on the 10% level. 
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To get an estimate of the quark mass ratio dependence, we have used two 
different sets of quark masses, one inspired by the chiral perturbation theory 



and the other by the quark masses extracted from the lattice [34( . The two sets 
have the following ratios 

1:1.83:36.7 (Mi) 

1 : 2.32 : 45.0 (M 2 ) ' 1 ' 

For each set we perform 5 simulations with ever smaller absolute masses, see 
Fig. [ini This data is fitted to ([521 to extract the chiral condensate. The results 
for the two sets agree on the la level, and we can argue that the chiral condensate 
depends only weakly on the quark mass ratios, given that the latter vary by 
roughly 25%, see ([54]) . This is as it should be since the exact chiral condensate 
does not depend on the quark masses at all. From an operational point of view, 
the robustness against quark mass ratios^ makes the chiral condensate a good 
quantity to set A. 

Solving (1531) we find that the lambda parameter is given by 

_ f 401(5)(40)(15)MeV . . 

1 ~ { 389(6)(40)(15)MeV ' 

where the errors follow from the fit, (qq)o v and the systematic on %. This leads 
to an overall error of 44MeV, or roughly 11%, and is strongly dominated by 
the one-loop result for (qq)^. Since we run at one loop in the self-consistency 
equation (l53l) . such a large error is certainly realistic, if not underestimated. 

We found that running at two-loorF^I gives results consistent with (f55|) . Usin 
the prescription of |9j , the scale and the mean instanton size for the IILM 



sing 



/ 598(65) MeV f 0.33(3) fm 

~~ \ 580(64) MeV ' pAi ~ \ 0.34(4) fm ' [bb > 

This is in very good agreement with the precision study [9(. Given that both 
works use chiral properties for the calibrations, the nice overlap is probably not 
totally unexpected. 

Note that (|55l) is a prediction for the lambda parameter with 3 active quark 
flavours. To compare our result with experimental data we run down the cou- 
pling constant a^f s = 0.117(2) [l6[ from Mz to and convert it to a lambda 
parameter. This is a rather big difference in scales and it is appropriate to use 
two-loop running, although not entirely consistent when we compare it to the 



18 i.e. taking the limit from different directions in quark mass space. 

19 We use /3-functions and anomalous dimensions from the MS scheme, [49| . since we do 
not know them for PV regularisation. However, both regularisations are thought to give 
roughly similar results, for instance (qq)^ v and (qq)^ 5 agree on the 3% level at one-loop and 
/i = 2 GeV, and for the purpose of estimating errors in the one-loop running this procedure 
should be fine. 

20 Remember that in the unquenched case the instanton size is fairly independent of the 
quark masses. 
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0.0034 




Figure 16: Computing the topological susceptibility allows for the extraction of the chiral 
condensate by using chiral perturbation theory, 1521 1, The rationale is the same as used in 
recent lattice studies to extract the chiral condensate 0, @, 01 • In order to get a rough estimate 
on the systematic error introduced by the chiral limit, two sets of masses have been used. The 
upper plot corresponds to the set Mi and the lower plot to M2, as given in 1541 1. The chiral 
condensate for both mass ratios agrees on the la level, and we conclude that it depends only 
weakly on the quark mass ratios. 
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one-loop result (|53T) . To deal with threshold effects, we use the Mathematica 
package RunDec, [4j| . The conversion between the MS and PV lambda parame- 
ters is given by [26| . Q 



Apv = AMsCxp( v 22 _ 4jV//3 j • (57) 

This leads to = 325(40) and the IILM result agrees on the la level. Trusting 
the perturbative running down to the rather low scale /i\ is a leap of faith. 
However, earlier studies have seen good agreement between IILM and lattice 
predictions for physical quantities, such as meson masses, and so the agreement 
between the lambda parameters might not just be a fluke. 

To determine the physical quark masses, we will use (|52j) rewritten in terms 
of the pion mass 

2 , 2 m u m d /1\ J (77.4 MeV) 4 



{m u + m d f \mj \ (75.9 MeV) 

We used = 135 MeV and /„ = 93 MeV. Together with the fits, Fig. [JJ 
we can compute the corresponding quark masses. We convert them into MS 
masses at 2 GeV, run at one-loop, in order to compare them more easily with 
other sources. Our results are 



mr^ = 0.&GeV) = \ ^ ™?A S (MeV) , (59) 



2.2(2) 4.0(4) 80(9) 
1.9(2) 4.4(4) 87(10) 

m P(^ 2 GeV, = {; : « «W -("'(MeV,. (00) 

The errors include an estimate from the 2-loop running. These masses compare 
well with the particle data group masses [l5j], i.e. m u = 1.5 — 3.3 MeV, ijid = 



3.5 — 6.0 MeV and m s = 70 — 130 MeV, and to the lattice masses [34|, i.e. 
m u = 1.9(2) MeV, m d = 4.4(3) MeV and m s = 87(6) MeV. 

Very large volume simulations are expensive even in the IILM. We have seen 
in section [5] that the instanton density becomes independent of quark masses in 
the chiral limit. Therefore, in the physical region of parameter space that we 
are considering the volume is directly proportional to the number of instantons 
in the box. The computation of the fermionic determinant is the most costly 
part of the simulations. In a naive implementation it would scale as 0(V 3 ). 
However, in section l4~3l we were able to drastically reduce this cost to 0(V 2 ) 
by implementing fast update algorithms for the Cholesky decomposition. The 
absolute scale for the volumes is set by considerations regarding hnite size effects. 
To have these under control the lightest propagating particle, in our case the 
pion, should fit into the box. For masses beyond the chiral limit the instanton 
density increases, see again section [5] However, the pion mass increases too and 
it turns out that for the same level of control over finite size effects the system 
size can be smaller, i.e. larger quark masses are computationally cheaper. 

We have studied the thermodynamic limit on four volumes, in the range 2 < 
Lm, < 3. Even though the data has displayed a nice scaling with the volume, it 
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2000 4000 6000 8000 10000 12000 14000 16000 

V (A" 4 ) 



Figure 17: The simulations are performed for different simulation boxes. The average of the 
smallest Dirac eigenvalue, (A m ; n ), is smaller than the quark masses for all but the smallest 
simulation box. 



is important to check whether the thermodynamic limit was consistent. To this 
end we will run large volume simulations, Lm, € [2.11,3.7], for the particular 
set of physical masses inspired by chiral perturbation theory: in dimensionless 
units, m u = 0.00546, m d = 0.01001 and m s = 0.2002. This will allow us to 
estimate the systematic error introduced by performing the thermodynamic on 
the set of smaller simulation boxes. 

The thermodynamic limit on the topological susceptibility turns out to be 
rather insensitive, see Fig. 1181 However, the mean instanton size does not con- 
verge to a constant even for the largest volumes, see Fig. [TH1 It is a rather 
lucky fact that the mean instanton size does not vary much in absolute terms. 
The slope is clearly decreasing and we might estimate the convergence to oc- 
cur somewhere in the range p £ [0.68,0.66]. A fit to p — p^ + aV~ 25 gives 
Poo = 0.6720(5) A -1 = 0.33(3) fm, in good agreement with the phenomenologi- 
cal value. The instanton density turns out to be n = 1-7(7) fm _ , and like the 
topological susceptibility displays a nice thermodynamic limit. 



7. Conclusion 



With the discovery of the new non-trivial holonomy calorons 2?1 3(J, 31 1 



there is renewed interest in studying the role of non-trivial field configurations in 
QCD, especially their role in the confinement/deconfinement phase transition. 
A lattice based approach [2l[ is well suited for the pure gauge sector because 
it easily incorporates many-body instanton interactions in the classical action. 
However, the introduction of fermions will be plagued by the same problems that 
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Figure 18: The fluctuations of the topological charge (Q 2 ) show a nice linear dependence with 
the volume box V, as it should be for an extensive quantity. Applying the thermodynamic 
limit to the 4 smallest volumes yields a topological susceptibility that agrees on the la level 
with the corresponding result using all available volumes. The mean instanton number ranges 
from (N) rj 200 to (AT) »j 1600. 
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Figure 19: Even for the largest volumes the mean instanton size p is still decreasing and does 
not seem to converge to a constant. We are rather lucky that, although the effect is clearly 
systematic, the variation of p is small in absolute terms. 
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full lattice computations face. Most notably, the computations will become very 
costly. 

A different approach, pioneered by Shuryak, Diakonov and Petrov, is to 
formulate the instanton liquid in the continuum, see for instance [39]. This 
approach suffers from the fact that it is not straightforward to include many- 
body interactions. Incidentally, only two-body effects are taken into account. 
The strongly localised profiles of instantons and the, a posteriori, fact that the 
instanton ensemble is rather dilute make this a viable working premise. The 
big advantage of the continuum formulation is the ease with which quarks can 



be incorporated, see also [5l|. From this perspective, both the lattice and the 



continuum models complement each other rather well. Meanwhile, different 
publications have investigated the confining nature of other backgrounds, such 
as regular gauge instantons and merons 32, 3(|, [H(J at zero temperature. 

So far the continuum models used explicit analytic formulas for the inter- 
actions. They have been obtained through asymptotic considerations and fits 
to numerical evaluations of the classical action. We noted that these formulas 
do not possess a thermodynamic limit at finite temperature. More importantly 
perhaps, the more complex moduli-space of the non-trivial holonomy calorons 
probably demand a more systematic approach. In this paper we have set up 
a framework which we believe is numerically well-defined, can be extended to 
more complicated backgrounds and does not suffer from the parametrisation 
bias introduced implicitly through analytical formulas motivated by symmetry 
arguments and fits. The price to pay is a larger numerical overhead because 
evaluation of the interactions through look-up tables and asymptotic matching 
formulas is computationally more expensive than through simple fitting formu- 
las. 

We have found that the analytic formulas of |39j agree very well with our 
interactions at zero temperature. Especially for the case of equal instanton sizes, 
where strong symmetry arguments support the analytic formulas, the agreement 
can be seen as a validation of the numerics. In general, however, the interactions 
of both schemes differ; the differences are especially pronounced for the quark 
overlaps because in this paper we use the full ratio ansatz whereas a sum ansatz 
is used in 39]. Shifting the point of view, we considered the formulas of [3!| 
to be another valid scheme; together with the streamline ansatz we studied the 
dependence of bulk properties of the IILM on the choice of these three rather 
different interactions and found that this introduces a systematic effect which 
depends on the quantity under consideration, but was generally rather large, up 
to 20%. 

The IILM has been shown to be compatible with the chiral properties of 
QCD, see for instance A key chiral property, the topological susceptibility, 
has not been studied extensively within the IILM, see however [46}. One reason 
might be that the IILM was so far set up for simulations in the canonical ensem- 
ble whereas the topological susceptibility is most naturally studied in the grand 
canonical ensemble. We have enlarged the Monte Carlo moves to incorporate 
insertion/deletion steps in order to simulate an open ensemble. Apart from 
technical problems related to book-keeping issues this is rather straightforward. 
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A major incentive for this work was to investigate the regime of physical 
quark masses. In order to deal with such light quarks, rather large volumes 
need to be considered. We dealt with this issue by, first, reducing the complex- 
ity of the algorithm from 0(N 3 ) to 0(N 2 ); this is achieved by rewriting the 
updates in a form suitable for fast matrix modifications. Secondly, we study the 
thermodynamic limit and monitor some bulk quantities to guarantee a consis- 
tent large volume extrapolation. 

The topological susceptibility is easy to compute in the IILM and has been 
studied extensively on the lattice. It represents a natural candidate to hx units 
in the IILM. For quenched simulations we found rather good agreement be- 
tween the IILM and the lattice in this way. We are, however, mainly interested 
in the unquenched case. Instead of using the topological susceptibility directly, 
we decided to use the chiral condensate, extracted through the topological sus- 
ceptibility and chiral perturbation theory, to set units. The reason we decided 
against a direct use of the topological susceptibility is that it depends strongly 
on the quark masses. We found that the chiral condensate has a very weak de- 
pendence on the chiral limit, and use it to set the units in the unquenched sector. 
We achieve good agreement with previous work and also with experimental data 
on the strong coupling constant a s . Using chiral perturbation theory, we are 
able to determine physical quark masses. These turn out to compare well with 
experimental bounds and lattice simulations. 

Finally, we investigated the uncertainties introduced by the large volume 
extrapolations, and found that our procedure, of bounding the volume in such 
a way that the quark masses are smaller than the smallest Dirac eigenvalues, 
allows for a systematic thermodynamic limit. 

In a further publication we will use these input parameters to study the 
IILM at finite temperature. 



Acknowledgements 

We are very grateful for many informative discussions with P. Faccioli and 
E.P.S. Shellard. Simulations were performed on the COSMOS supercomputer 
(an Altix 4700) which is funded by STFC, HEFCE and SGI. This work was 
supported by STFC grant PPA/S/S2004/03793 and an Isaac Newton Trust 
European Research Studentship. 



Appendix A. Gluonic Interactions 

The ratio ansatz for an instanton-anti-instanton pair is defined by 

A a = V^dMx, {X UP1 }) + O a % v d v Il 2 (x, {X 2 ,p 2 }) 

" l + ni(a;,{a;i,pi})+n2(s,{z2,P2}) ' 

where H(x, {y,p}) = ^ and O = 0\0 2l with 0% the respective colour embed- 
dings. A global colour rotation has been performed to bring the gauge potential 
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into this form, which is irrelevant since the action is gauge invariant. Instanton- 
instanton and anti-instanton-anti-instanton pairs differ by having either only fj 
or 77 in the above formula. A brute force computation then gives 



= I + ( Tr0t ° + inor,)^)J + (vOv) P » pv i^ 

+ {T]O t 0^ wva J wva + {r\Ori) a ^ ap {riOr\) lilJ f icr K il(nJ<J . (A.2) 
The different terms have the following form 

1 = (1 + ni 4 +ri2)2 [(d.AniX^n!) + (0 M a,n 2 )(0„a,n 2 )] 

- (1 + rii 8 +n2)3 [(^a,n 1 )(^n 1 )(s,n 2 ) + (^n^n^no 

+2(^^n 1 )(a /1 n 1 )(5 i ,n 1 ) + 2(^^n 2 )(^n 2 )(^n 2 )] 

+ (1 + ni 4 + n2)4 [3(a„ni^ni)(a„ni^n 2 ) + 3(^n 2 a M n 2 )(^n 2 a Al n 1 ) 

+3(^n 1 a M n 1 ) 2 + 3(d^n 2 d^n 2 ) 2 + 2(^n 1 ^n 1 )(^n 2 ^n 2 ) 

+(a M n 1 a AJ n 2 ) 2 ] . (a.3) 



3 = (i + n 1 + n 2 )4 ( d ^W( d ^W • (A- 4 ) 



^ = (i + n 1 4 +n a )» (0 ^ ni)(0 ^ na) 

+ (1 + ni 4 + n2)3 [(^^no^na^na) + (a^n 2 )(a ff n 1 a ff n 1 ) 

-2(9 Al a ff ni)(a iy n 2 )(9 CT n 2 ) - 2(d t ,u 1 )(d rT u 1 )(d,d a u 2 ) 
-2(9 Al a ff n 1 )(9 ff n 1 )(a,n 2 ) - 2(^n 1 )(d 1 ^n 2 )(d CT n 2 )] 

+ (1 + n 4 + n2)4 [-(^n 1 )(^n 1 )(9 ff n 2 a CT n 2 ) - (^n^n^n^n!) 

+3(d^u 1 )(d lJ u 2 )(d a ii 1 d a u 1 ) + 3(^ni)(^n 2 )(9 ff n 2 a CT n 2 ) 

+3(d ll u 1 )(d v ii 2 )(d <r u 1 d a u 2 )] . (A.5) 



+ (1 + rii 8 +n2)3 [(^n 2 )(a p a,n 1 )(a (T n 2 ) + (^n 1 )(a p a,n 2 )(a ff n 1 )] 

+ (1 + ni 8 + n 2 )4 (^ni)^)^^)^) . (A.6) 
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■W = (1 + Ui + n2)4 (^ n i)(^ n 2)(^n 1 )(^n 2 ) . (A.7) 



^ CT = (i + it + n 2 )4 (^n 1 )(^n 2 )(^n 1 )(a g n 2 ) . (A.8) 

Appendix A.l. Exact Interactions 

When computing the look-up tables, we use global translations and rota- 
tions in R 4 to place one instanton at the origin and the partner &t y' 4 = R = 
yfR^R^ = \y !l — y l2 \, where y % are the instanton centres. The rotation will 
reemerge in contractions of R^ with the colour structure, as we will now see. 
The relation between the position vector i? M and R'^ = (0, 0,0, R) is given by 
the following rotation matrix 

Rfi = Ofu/Rv > 

<V = ^, (A.9) 

and the other components of the rotation matrix are irrelevant. 

Note that, with the choice of R' the integrands are 0(3) symmetric in the 
subspace orthogonal to the 4-direction. Denoting the arguments of the 't Hooft 
potentials U(x,{y,p}) by x^ and x^ = x^ — we can extract extract the 
R^ from the integrands with help of the following formulas, which we order 
according to the tensor structure of the ^-dependence on the integrand. 

J x lt = O l * J x' 4 . (A. 10) 

j x^x v = <J M „ J x? + 0^0„i J {xf - x?) . (A.11) 

j x^x v x K = (S^O k4: + 5 KfJt O vi + S vk O^) J xfx' 4 (A.12) 
+ O^O vi O Ki J(x'i - 3x'?x' 4 ) . (A.13) 

J x^x v x K x s = (S^6 kS + 5p K 5 v s + S^sSkv) J x?x$ ( A - 14 ) 
+ [S^O^Osa + porm.) J (xfxf - xfx'f) (A.15) 
+ O^O uA O kA Osa j (x' 4 - 6xfxf + Sxfx'l) . (A.16) 
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Terms with x can be constructed from these. Incidentally, splitting the 
different integrands according to the above formulas is the most stable procedure 
numerically. Taking into account the antisymmetry of the 't Hooft symbols, we 
end up with the following integrands. 



/ = 



+ 



(i + n! + n 2 ) 2 

8 

(i + n 1 + n 2 )3 

4 

(l + IT+H,) 4 



[(ni') 2 + 3(ni/r) 2 + (i^') 2 + 3(i^/f) 2 



XX 



2W{(W 1 y + 2n 2 '(n 2 ) 2 + (n'/nin 2 + n;n 2 'n 2 ) 



, XX ^ 



i tt' ^2 



i2(ii' 1 r + i2(n 2 ) 4 + 8(n[y(ii' 2 y + 4(-;rrin 2 

rr 



+i2(n' 1 ) 2 (^n' 1 n 2 ) + i 2 (^n;n 2 )(n 2 )^ 



. (A.17) 



Note that to achieve good numerical precision, we need to subtract the 
one-instanton integrands from the above before performing the numerical inte- 
gration. 



J - (i + n 1 + n 2)4 ^ 2 ^ 2 - (A.i8) 



T —51+ R H R v j (A lq \ 



/w (l + ni + n 2 ) 2 



UU'I - (U\/r))(U' 2 /r) + ^Wr)(n 2 ' - (Il' 2 /f)) 



+(n;/r)(n 2 /r) + X X ^(n';n 2 ' - n'/(n 2 /r) - (n;/r)n 2 ' + (u[/r)(W 2 /f)) 

rr rr 



+ 



1 



(i + Hi + n 2 ) 3 



[4((n' 1 /r)(n 2 ) 2 + (n' 1 ) 2 (n' 2 /f)) 



+ 



+^(4(n'/ - (ni/r))^) 2 - siii^iu'jf)) 
+||(4(n' 1 ) 2 (n 2 ' - (n 2 /r)) - scni/r)^) 2 ) 
+^(-8^((n'/ - (n[/r))(n' 2 y + (n;) 2 ^' - (n' 2 /f))) 

rr rr 



-sn'/nin, - snin^'n',] 



i 



(l + ni + n,) 4 



- 4 ^(n' 1 ) 2 (n 2 ) 2 -4^(n' 1 ) 2 (n 2 ) 2 



„'2 



+i2^n' 1 n 2 ((n' 1 ) 2 + (n 2 ) 2 + ^n;n 2 ) 

rr " rr 



. (A.20) 
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(i + n! + n 2 ) 2 



+ K y x? (K/rm (n' 2 /f)) 
+ <(<- R )- x ? ™ {u >;iL» n'/(n 2 /f) (ni/r)n 2 ' + (n;/r)(n 2 /r)) 



rr 



+ 



1 



(i + ni + n 2 



.T-^ 



(4(n? - (n;/r))(n 2 ) 2 - 8(n;) 2 (n' 2 /f)) 



(4(n' 1 ) 2 (n 2 '-(n 2 /f))-8(n' 1 /r)(n 2 ) 2 ) 



+ <K fl) - (-8^((nr - (n;/r))(n 2 ) 2 + (ni) 2 (n 2 ' - (n' 2 /f))) 

rr rr 



+ 



(i + n! + n 2 ) 4 



-8n'/n;n 2 - 8n;n 2 'n 2 ] 
-4 ^ "^ (ni) 2 ^) 2 - 4^— ^(n;) 2 ^,) 2 



, x 4( x 4 ^) ^'l n'n' 



+i2 -"-« 1" - 1 n^n.an;) 2 + (n' 2 ) 2 + ^n;n 2 ) 



rr 



rr 



(A.21) 



Urn./ „^ J /a 



i? 2 



(A.22) 



= (analytically) . 



(A.23) 



T~.l1 V 2 



(i + ni + n 2 ) 2 

- i?) 2 - xf 



^(ni'-(ni/r))(i^/f) 



+- 



x?R 2 



(i + n! + n 2 )3 



(ni/r)(n 2 ' - ( n 2 /r)) 7 i^(n' 1 ' - (n;/r))(n 2 ' - (n' 2 /f)) 

(x> 4 - Rf - xf 



f2 V"l/'/V" 2 // ( r ^2 



r /2 _ /2 



x' 2 i? 2 



■(nl/r)^) 2 



+^ 2 -((n'/ - (n;/r))(n 2 ) 2 + (n;) 2 (n 2 ' - (n 2 /f))) 

(rrj z 



+ 



(i + n! + n 2 ) 4 



x' 2 R 2 
(rf) 2 



(ni) 2 (n 2 ) 2 



(A.24) 



x ,2 R 2 



R R 2 

■W - <W p R2 a (1 + ni + n 2)4 {rf) 2 



(n;) 2 (n 2 ) 2 . 



(A.25) 
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K^pva — [(^i/'Jpcr + 5^p5 va + 5 pab~ pv)X-{ X 



12 12 
2 



RpRg { { , „>2 12 12 l2\ | r Rfj,Ri> , 12 12 12 I2\ 
' ft2 \\ X i U ' X l x l x 2 ) + Q pa ( X * X ~L ~~ X l X 2 ) 



+ (Sfip—^ h Sp a ^ 2 P + 5 up —jp h <W ^ 2 P ) (x' 4 (x' 4 - -R)^'] 2 — xfx'2) 
H ~ — ( X 4 i X 4 ^) 3*^1 ^2 (^4 R) x l 

-MK - *K 2 )] (1 + ni 2 +n2)4 (^(n;) 2 (n 2 ) 2 . (A.26) 

Appendix A. 2. Asymptotic Interactions 

As explained in the main text, the small separation asymptotic formulas 
get contributions which have the same functional form as those for the large 
separation asymptotics; the difference lies in the integration limit. We will 
therefore start with the large separation formulas and leave the integrals explicit. 

Appendix A. 2.1. Large Separation 

The upper integration limit z follows from variable substitution and has the 
the following for integration over Ii, with I 2 held fixed, 

z 2^i+ik r 2. (A27) 
Pi 

Apart from the dependence of z on II, the rational form of the 't Hooft 
potential allows for a complete factoring out of II under the above mentioned 
variable substitution. For the large separation formulas it is understood that 
z 1 — > oo because the initial integration variable r extends to infinity. 

The integral over / contains terms that do not mix the 't Hooft potential IT 
and II2 except for the denominators. At zeroth order in our expansion, these 
terms can be transformed to exactly match the single instanton contributions 
by exploiting scale invariance. Remembering that we actually subtract the 
one-instanton contributions to get the interactions, we can neglect these terms 
altogether. We then end up with the following formulas. 



f r ™ 2 2 dnUdnU [* S 5 d S 

J I = 1 ^ p t+wl FTi7 + sym ' (A - 28) 
r 2 2 dpUdpU f z s 5 ds 

J J = 167Tp OTWJ WTW + sym - ( 9) 
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f T 2 2 df,d u n f 

J^ = 167Tp (YTWJ 



(i + n) 2 7 (s 2 + i) 3 
/ 2 2A «n xs 2 2 (W)\ r s^ds 

At zeroth order, partial integration and the antisymmetry of the 't Hooft 
symbols can be used to simplify 

w ~* (i + ni 8 + n 2 ) 4 (^ni)(^n 2 )(^n 1 )(^n 2 ) , (A.si) 

with asymptotic behaviour 

/ W = ^PX ^l^ f + ^ • (A.32) 



/* T ,2 2c (d„m{d a Tl) f z s 5 ds 

J J ww = ^ l [ 1 ^3 / (^TTF + sym - (A - 33) 

For K wva no 't Hooft symbols can be used to exchange the index pairs 
f) <H> (p, a), and so we cannot simplify with a symmetry argument anymore. 



Ik 4^2. (d p u 2 )(dM p a 5 ds 
J wv ° ~ n Pl °^ (i + n 2 ) 3 J (s 2 + i) 4 

/ 



+ 47TM- (1 + ni)3 y ^ 

Appendix A. 2. 2. Small Separation 

As explained in the main text, the small separation asymptotic formulas get 
contributions from the large asymptotics. Also, in this case we have performed 
a global translation so that the instantons sit at ±i? M /2. Therefore, in the large 
separation formulas we need to put z 2 = i^(i?/2) 2 . 

We now turn to the proper small separation asymptotic formulas that encode 
the repulsion through the gauge singularity. We will again introduce an explicit 
upper limit for the integrals; abusing notation we will use the same letter as 
before, but here the meaning becomes 

R 2 R 2 

z2 = pTT7 2 > zf = 7C (A - 35) 



To derive these formulas, we approximate the arguments x^ ± i? A ,/2 



x 



We have, therefore, explicitly restored 0(4) symmetry, which can be exploited 
to set several integrals to zero. Eventually, we arrive at 
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I = 384tH 



4 i 4 

PI + Pi 



ds 



2 2 

PIP2 



(p? + p!) 2 (p? + p!)V y,*(* 2 + i) 3 
pf + pi + + pfpl + p\pI f ds 

(pi+ptr y*s(* 2 +i) 4 



s(s 2 + l) 4 7 22 s(s 2 + l) 4 



(A.36) 



J = 64tt 



4 4 
2 P1P2 



ds 



{pI + pIY h ^ 2 + i) 4 



(A.37) 



96tt 



2 P1P2 f d s 



(Pl+Pl) 2 J z s( S i + 1) 



2 PiP2+3pfp|+3p^ 



1927T 



2^6 



2 p?p! 



ds 



{pl + pl) 2 J z s(s 2 + l) 3 



-32tt 



(Pi + pl) 



2\4 



ds 



s(s 2 + l) 4 



(A.38) 



-32tt 



2 P1P2 f ds 



{Pi + Pl? h s(s 2 + l) 



-32tt 



2 2 

2 PIP2 



ds 



(p 2 i + p 2 2) 2 J z S ( s 2 + iy 



(A.39) 



J p pis a • 



(A.40) 



K, 



4 4 
P1P2 



ds 



(A.41) 



Appendix B. Fermionic Interactions 

The Dirac overlap matrix elements are given by 



Tia = I d 



(B.l) 



Note that iTr({7rj") = iup is the colour four- vector used for instance in |39j . 
After some straightforward algebra, we find that Ip has the following form 



1 



(1 + u z + u A )(i + n 7 )3/2(i + n A )3/ 2 



l + Ei 



{dpUjdpU^dpUA + (d^UAd^dpU! . (b.2) 
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Appendix B.l. Exact Interactions 

Using the same rotations (|A.9[) as for the gluonic interactions to marry the 
space-time with the colour indices, we get 



j = Rfi =3 

R (i + n, + u A )(i + n/)3/ 2 (i + n A ) 3 /2 

{^(n^ + ^W) 2 ^^}. (B.3) 

Appendix B.2. Asymptotic Interactions 
Appendix B.2.1. Large Separation 

In order to get rather simple formulas, we make the following additional 
simplification 

1 + n 7 + II 4 \ ] + ^ : ^egration over H, 

\ 1 + II A : Integration over Il A v ' 

Given these further assumption, we can proceed as for the gluonic interac- 
tions. Finally, caution needs to be taken in the case of an anti-instanton because 
it sits at — i? M so that O^Ha generates an extra minus sign. 



2 2 u A d n A 

07T pj- 



(i +1X4) 3 / 2 y (s 2 + i) 7 / 2 

9 , d B Ui f ZA s 4 ds , x 

Appendix B.2.2. Small Separation 

At zeroth order, i.e. x ll ±R fl /2 —} x^, the contribution to Ip vanishes because 
of 0(4) symmetry. It turns out that the large separation asymptotics falls off 
too quickly as R — > 0. However, this is not important because in this regime 
the gluonic interaction is dominant. 



Appendix C. Cholesky decomposition update 

In this appendix we look in detail at how the structure suitable for the 
Cholesky decomposition update comes about. We will also see that insertion 
can be performed faster whereas deletions will be the most costly. 
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Appendix C.l. Canonical Moves 

Upon updating instanton /, we have that T — > T+AT, with (AT)y = <5ij£j. 
This induces the following changes 

(T^T)ij -+ (T^T) ij+ T^*+^T Ij+ ^*, (C.l) 
V>« = T*, (C.2) 

4>i = Ci + V'i, (C3) 



+ ^u^k^lj + TikZkdij + SuSji , (C.5) 

V>i = j^ (T ^' (c - 6) 

& = tf/i^+Vi, (C7) 
T t T _> T^T + 4>(^ - ^ . (C.8) 

Changes in an anti-instanton will have analogous formulas. 

Appendix C.2. Insertion 

Wc focus on inserting an instanton. Insertion of an anti-instanton is then 
similar. Since in the code we always add an instanton at the end of the arrays, 
an insertion corresponds to adding a row to T. 



T ( ^ ) . (C.!)) 



^ - (?)( Tt ° = (??t ?!)■ (c - l0) 

On the level of the Cholesky decomposition 

L 

X 1 



L -+ = ( f J ) , (C.ll) 



* . -(» J), (C, 2 ) 

Remembering that the insertion also adds a mass term in the diagonal, we 
have to solve the following system 

r LD X = n (cu) 
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which can be solved in 0(N 2 ) by using backsubstitution. The case for T^T is 
simply given by 

T t r ^( T t C ) f ^ ) =^ t T + ^ t , (C.15) 
which is a rank-1 update. 

Appendix C.3. Deletion 

We focus again on an instanton. Deletion will be a two step process. We 
first delete the last instanton and then swap it with that instanton that has been 
selected for deletion. The swapping is similar to a canonical move, where now £ 
is given by the difference between the last instanton and the selected instanton. 

The proper deletion part is given by 

TT t _> 
L -> 
D -> 

The T^T part is again simply related to a rank-1 update because, upon 
rearranging the result from the insertion part, we get 

T^T -> T^T — . (C.19) 



TT 1 


L 
1 

D 




(C.16) 
(C.17) 
(C.18) 



Appendix D. MS to PV 

Operators with anomalous dimensions run, and for mass independent renor- 
malisation prescriptions they depend on the scheme. The IILM makes predic- 
tions within a subtraction scheme that uses Pauli-Villars regularisation. How- 
ever, the lattice results have been quoted in MS, and so we need to work out 
the relation between the two. 

It is not hard to convince ourselves that the quark masses run inversely 
to the chiral condensate: note that chiral perturbation theory, (|58jl. relates 
the topological susceptibility to the pion mass and decay constant, which are 
physical quantities; it also relates the chiral condensate and the quark masses 
to the topological susceptibility through and therefore the renormalisation 
scheme dependence must exactly cancel among the two. 

Eventually, we will also relate the quark masses of the two schemes, and so 
here we will focus on mass renormalisation. We will only work at one-loop^], 
i.e. we will need to evaluate (Fig. ID.20[) in both schemes. 



21 Maintaining manifest gauge-invariance in Yang-Mills theories using Pauli-Villars regular- 
isation is not straightforward beyond one-loop. 
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Figure D.20: Feynman diagram needed to compute the difference between the MS and PV 
scheme at one-loop. 

This is a textbook computation, [38| . After subtracting off the divergences, 
we end up with 

Spv = gc(3) {-2™ + \i> + J M 2m - (1 - X )i>) ^ xm2 _f {1 _ x)p2 } , 

(D.l) 

and 

^ = ^C(3) { -m + \i + / ^ - (I - X )i>) In xm2 J {l „ x)p2 } • 

(D.2) 

Relating both through the pole mass and using that, in our case C(3) = 4/3, 
we get 

to ms = mpv(l- r^r). (D.3) 
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